|
|
\[\textbf{Lineu Alberto Cavazani de Freitas}\] \[\textbf{Prof. Cesar Augusto Taconeli}\] \[\textbf{Prof. José Luiz Padilha da Silva}\]
Análise Comportamental de Ovelhas Submetidas a Intervenção Humana usando GAMLSS
Material Rbras: Número de Mudanças de Postura de Orelha e Cauda
O conjunto de dados é proveniente de um experimento sobre o comportamento de ovelhas, conduzido na fazenda experimental INRA La Fage, Roqueford, França, em setembro de 2015 com o objetivo de verificar o efeito de linhagem genética, escovação e isolamento nas respostas comportamentais dos animais (TAMIOSO et al., 2017). Na ocasião, vinte ovelhas classificadas como reativas ou não reativas ao isolamento social temporário foram submetidas à escovação por um humano familiar. As ovelhas tinham 15 meses de idade, eram não gestantes e não amamentavam quando foram observadas.
|
|
O experimento foi conduzido em três sessões experimentais: na primeira tinha-se uma grade de metal separando o animal testado dos demais animais, sem distância entre eles. Na segunda havia duas grades de metal separando os animais a uma distância de 1,7 metros, ou seja, foi imposta a condição de isolamento social. E na terceira sessão, os animais voltaram a ser separados por apenas uma grade:
|
|
As sessões de testes ocorreram dois dias após a fase de adaptação dos animais ao equipamento e aos humanos e, em cada sessão, as ovelhas foram observadas em 3 momentos distintos: fase de pré escovação, com duração de 2 minutos e 30 segundos; fase de escovação, com duração de 3 minutos; e pós escovação, com duração de 2 minutos e 30 segundos. O delineamento para um animal pode ser represntado da seguinte forma:
|
|
Os dados coletados dizem respeito ao número de mudanças de postura dos animais e à proporção do tempo em que os animais permaneceram em determinadas posturas, tratando-se então, de um conjunto de dados com múltiplas respostas em que não há observações independentes, já que cada animal contribui com nove medidas. Portanto, há a necessidade de incorporar as correlações entre as medidas num mesmo animal e do animal dentro de cada sessão experimental, além da correlação entre as respostas.
Foram avaliados os efeitos de:
Sessão: Fator de 3 níveis que indica a sessão experimental (Se1, Se2, Se3).
Momento: Fator de 3 níveis em que indica o momento experimental (antes, durante, depois).
Linhagem: Fator de 2 níveis que classifica os animais como reativos ou não reativos ao isolamento social temporário.
setwd("~/Área de Trabalho/PESQUISA/IC_nova/1. Dados")
## Pacotes
library(ggplot2)
library(gridExtra)
library(gamlss)
library(knitr)
######################################################################
## Dados
count <- read.csv2("DataCont.csv", header = T, sep = ',')
## Alterando os contrastes default para os contrastes de interesse
count$sessaoR <- count$sessao
count$sessaoR <- relevel(count$sessaoR, ref = 'Se3')
contrasts(count$sessaoR) <- 'contr.helmert'
contrasts(count$sessaoR)[,1] <- contrasts(count$sessaoR)[,1]/2
contrasts(count$sessaoR)[,2] <- contrasts(count$sessaoR)[,2]/3
## Criando um fator combinando animal e sessao
count$animals <- factor(count$animal):count$sessao Para o ajuste dos modelos de regressão foram consideradas as distribuições:
Para cada uma das 7 distribuições ajustou-se o modelo no qual era considerado sessão, momento, linhagem, as interações duas a duas e efeitos aleatórios a nível de animal e animal dentro de sessão para a média.
Nas distribuições com parâmetros referentes à inflação em 0 também foram acrescentadas as mesmas covariáveis, com exceção dos efeitos aleatórios.
Todas análises foram feitas utilizando o pacote gamlss do software estatístico R.
modelo_norelha <- function(familia){
gamlss(norelha ~ sessaoR + tempo + linhagem
+ (sessaoR + tempo + linhagem)^2 +
re(random = list(animal=~1, animals=~1), method = 'REML'),
nu.formula = norelha ~ sessaoR + tempo +
linhagem+ (sessaoR + tempo + linhagem)^2,
data = count, method = RS(50),
family = familia)}
mPO <- modelo_norelha(familia = PO())
mNBI <- modelo_norelha(familia = NBI())
mZIP <- modelo_norelha(familia = ZIP())
mZINBI <- modelo_norelha(familia = ZINBI())
mZAP <- modelo_norelha(familia = ZAP())
mZANBI <- modelo_norelha(familia = ZANBI())A escolha da distribuição utilizada baseou-se nos valores observados para as medidas AIC, BIC, Deviance, Verossimilhança, graus de liberdade e por fim, verificou-se se houve convergência.
| Modelo | AIC | BIC | Deviance | logLik | df | Conv | Iter |
|---|---|---|---|---|---|---|---|
| PO | 1191.028 | 1344.848 | 1094.6787 | -547.3393 | 48.17479 | TRUE | 5 |
| NBI | 1054.742 | 1206.695 | 959.5615 | -479.7808 | 47.59021 | TRUE | 7 |
| ZIP | 1127.398 | 1315.057 | 1009.8519 | -504.9260 | 58.77281 | TRUE | 17 |
| ZINBI | 2488.164 | 2687.903 | 2363.0517 | -1181.5259 | 62.55619 | FALSE | 50 |
| ZAP | 1137.740 | 1324.007 | 1021.0663 | -510.5332 | 58.33675 | TRUE | 12 |
| ZANBI | 1046.170 | 1247.390 | 920.1308 | -460.0654 | 63.01976 | TRUE | 10 |
A distribuição selecionada foi a Binomial Negativa Zero Ajustada.
Considerando o modelo ZANBI, foram realizados 3 reajustes para o modelo original:
Todos os modelos reajustados são modelos encaixados ao modelo original, sendo assim utilizou-se o teste de razão de verossimilhanças para selecionar o menor modelo possível que não difere estatísticamente do original.
mZANBIr1 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
(sessaoR + tempo + linhagem)^2 +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
nu.formula = norelha ~ sessaoR + tempo + linhagem,
family = ZANBI(), data = count)mZANBIr2 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
nu.formula = norelha ~ sessaoR + tempo + linhagem,
family = ZANBI(), data = count)mZANBIr3 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
family = ZANBI(), data = count)## Modelo AIC BIC Deviance logLik df Conv P-val
## 1 Completo 1046.170 1247.390 920.1308 -460.0654 63.01976 TRUE -
## 2 Reajuste 1 1036.359 1212.034 926.3192 -463.1596 55.01976 TRUE 0.63
## 3 Reajuste 2 1022.671 1174.258 927.7194 -463.8597 47.47563 TRUE 0.95
## 4 Reajuste 3 1037.248 1172.871 952.2970 -476.1485 42.47563 TRUE 0.05
Conclusões
O modelo sem interações para nu não difere do modelo completo.
O modelo sem interações para mu e nu não difere do modelo completo.
O modelo sem interações para mu e nenhuma covariável para nu difere do modelo com covariáveis para nu.
Sendo assim, as covariáveis para mu e nu devem ser mantidas, porém sem as interações.
Para o número de mudanças de postura de cauda foi experimentado um modelo com os efeitos fixos de sessão, momento e linhagem na média e na dispersão. Além disso considerou-se o efeito aleatório a nível de animal na média.
mNBI2_sigma <-
gamlss(ncauda ~ sessaoR + tempo + linhagem +
re(random = list(animal=~1)),
sigma.formula = ~ sessaoR +tempo+ linhagem,
data = count, family = NBII(), method = RS(100))Parte importante da etapa de ajuste de modelos de regressão é a análise de resíduos. A análise de resíduos em modelos da classe GAMLSS é baseada em resíduos quantílicos aleatorizados, estes resíduos são de interpretação mais simples visto que, sob um modelo bem ajustado, possuem distribuição Normal.
## ******************************************************************
## Family: c("ZANBI", "Zero Altered Negative binomial type I")
##
## Call: gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
## re(random = list(animal = ~1, animals = ~1), method = "REML"),
## nu.formula = norelha ~ sessaoR + tempo + linhagem,
## family = ZANBI(), data = count)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.49220 0.07710 32.324 < 2e-16 ***
## sessaoR1 0.42218 0.10556 3.999 0.000105 ***
## sessaoR2 -0.16399 0.09251 -1.773 0.078580 .
## tempoDepois -0.43246 0.09736 -4.442 1.86e-05 ***
## tempoDurante -1.05623 0.11387 -9.276 4.36e-16 ***
## linhagems+ 0.28393 0.08580 3.309 0.001205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9154 0.2149 -8.912 3.43e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: logit
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.1780 1.0714 -3.899 0.000153 ***
## sessaoR1 -2.6613 1.0911 -2.439 0.016049 *
## sessaoR2 0.4330 0.7312 0.592 0.554710
## tempoDepois 1.4892 1.1487 1.296 0.197069
## tempoDurante 2.7598 1.0831 2.548 0.011973 *
## linhagems+ -1.1191 0.6379 -1.754 0.081693 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 180
## Degrees of Freedom for the fit: 47.47563
## Residual Deg. of Freedom: 132.5244
## at cycle: 14
##
## Global Deviance: 927.7194
## AIC: 1022.671
## SBC: 1174.258
## ******************************************************************
## ******************************************************************
## Family: c("NBII", "Negative Binomial type II")
##
## Call: gamlss(formula = ncauda ~ sessaoR + tempo + linhagem +
## re(random = list(animal = ~1)), sigma.formula = ~sessaoR +
## tempo + linhagem, family = NBII(), data = count,
## method = RS(100))
##
## Fitting method: RS(100)
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.03866 0.14764 7.035 6.35e-11 ***
## sessaoR1 0.07220 0.17355 0.416 0.678
## sessaoR2 0.05521 0.14472 0.382 0.703
## tempoDepois -0.09011 0.21369 -0.422 0.674
## tempoDurante 1.66707 0.15827 10.533 < 2e-16 ***
## linhagems+ 0.16417 0.13724 1.196 0.233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.80294 0.38520 2.084 0.0388 *
## sessaoR1 -0.56679 0.42998 -1.318 0.1894
## sessaoR2 -0.07557 0.37539 -0.201 0.8407
## tempoDepois 0.55825 0.49029 1.139 0.2567
## tempoDurante 1.15304 0.42972 2.683 0.0081 **
## linhagems+ -0.00657 0.36109 -0.018 0.9855
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 180
## Degrees of Freedom for the fit: 27.95139
## Residual Deg. of Freedom: 152.0486
## at cycle: 13
##
## Global Deviance: 913.2412
## AIC: 969.1439
## SBC: 1058.392
## ******************************************************************
TAMIOSO, P. R.; BOISSY, A.; BOIVIN, X.; CHANDEZE, H.; ANDANSON, S.; DELVAL, E.; HAZARD, D.; TACONELI, C. A.; MOLENTO, C. F. M. Does emotional reactivity influence behavioral and cardiac responses of ewes submitted to brushing? Behavioural Processes, p. np, 2017.
|
|
|
|
|
|